In [3]:
function make_my_ode_closure(kappa, gamma, chi, epsilon)
function my_ode!(t, z, zdot, params=nothing)
zdot[1] = -.5kappa * z[1] + chi * z[3] * conj(z[2])
zdot[2] = -.5kappa * z[2] + chi * z[3] * conj(z[1])
zdot[3] = -.5gamma * z[3] - chi * z[1] * z[2] + epsilon
end
end
type my_odeparams
kappa::Float64
gamma::Float64
chi::Float64
epsilon::Float64
end
function eps_threshold(kappa, gamma, chi)
kappa * gamma /(4chi)
end
function my_ode2!(t, z, zdot, params::my_odeparams)
zdot[1] = -.5params.kappa * z[1] + params.chi * z[3] * conj(z[2])
zdot[2] = -.5params.kappa * z[2] + params.chi * z[3] * conj(z[1])
zdot[3] = -.5params.gamma * z[3] - params.chi * z[1] * z[2] + params.epsilon
end
Out[3]:
In [4]:
kappa = 10.
gamma = 100.
chi = 1.
et = eps_threshold(kappa, gamma, chi)
eps = 1.3 * et
Out[4]:
In [5]:
my_ode1! = make_my_ode_closure(kappa, gamma, chi, eps)
my_ode2_params = my_odeparams(kappa, gamma, chi, eps)
function my_ode3!(t, z, zdot, params=nothing)
zdot[1] = -.5 * 10. * z[1] + 1. * z[3] * conj(z[2])
zdot[2] = -.5 * 10. * z[2] + 1. * z[3] * conj(z[1])
zdot[3] = -.5 * 100. * z[3] - 1. * z[1] * z[2] + 325.
end
Out[5]:
In [39]:
tlist = linspace(0,10,2001)
z0 = zeros(Complex128, 3)
z0[1] += .1
hmax = .001
@time rk4solve(my_ode1!, z0, tlist, hmax, nothing)
@time rk4solve(my_ode2!, z0, tlist, hmax, my_ode2_params)
@time data_det = rk4solve(my_ode3!, z0, tlist, hmax, nothing)
Out[39]:
In [40]:
using Sundials
function wrap_complex(ode)
function real_ode(t, y, ydot)
ode(t, reinterpret(Complex128, y), reinterpret(Complex128,ydot))
end
end
@time Sundials.ode(wrap_complex(my_ode1!), reinterpret(Float64, z0), tlist)'
Out[40]:
In [41]:
function make_my_sde_closure(kappa, gamma, chi, epsilon)
function my_sde!(t, z, w, zdot, params=nothing)
zdot[1] = -.5kappa * z[1] + chi * z[3] * conj(z[2]) - sqrt(kappa)/2. * (w[1] + 1im*w[2])
zdot[2] = -.5kappa * z[2] + chi * z[3] * conj(z[1]) - sqrt(kappa)/2. * (w[3] + 1im*w[4])
zdot[3] = -.5gamma * z[3] - chi * z[1] * z[2] + epsilon - sqrt(gamma)/2. * (w[5] + 1im*w[6])
end
end
Out[41]:
In [53]:
tlist = linspace(0, 100, 20001)
my_sde! = make_my_sde_closure(kappa, gamma, chi, eps)
@time data_det = rk4solve(my_ode1!, z0, tlist, hmax, nothing)
@time data_stoch, w_stoch = rk4solve_stochastic(my_sde!, z0, tlist, hmax, 6, nothing)
Out[53]:
In [54]:
using PyPlot
In [56]:
plot(tlist, real(data_stoch[3,:]'))
plot(tlist, real(data_det[3,:]'))
Out[56]:
In [57]:
plot(tlist, real(data_stoch[1,:]'))
plot(tlist, real(data_det[1,:]'))
Out[57]:
In [58]:
plot(tlist, imag(data_stoch[1,:]'))
plot(tlist, imag(data_det[1,:]'))
Out[58]:
In [59]:
plot(real(data_stoch[1,:]'),imag(data_stoch[1,:]'))
Out[59]:
In [62]:
plot(tlist, real(data_stoch[1,:] .* data_stoch[2,:])')
plot(tlist, real(data_det[1,:] .* data_det[2,:])')
Out[62]:
In [ ]: